/*------------------<--  Start of Description -->--------------------\
| Power for two independent means;                                   |
|--------------------<--  End of Description -->---------------------|
|--------------------------------------------------------------------|
|--------------<--  Start of Files or Arguments Needed -->-----------|
| Arguments:                                                         |
|   - Required:                                                      |
|     std = pooled standard deviation                                |
|     n1,n2 = sample sizes for groups 1 and 2 respectively           |
|     y1 = true group 1 mean under null hypothesis                   |
|     min_y2 = smallest possible true group 2 mean                   |
|   - Optional:                                                      |
|     max_y2 = largest possible true mean for group 2                |
|     inc_y2 = increment value for range of y2                       |
|     alpha = type 1 error, e.g. .01 or .05, default=.05             |
|     sides = 1 or 2 for 1 or 2 sided test, default=2                |
|     plot  = 'P' for line printer plot of power vs y2               |
|             'G' for SAS/GRAPH plot of power vs y2                  |
|     unit  = units for mean, e.g. mg, lbs., cm, etc.                |
| Output: Power for true group 1 mean of y1 vs true group 2 means    |
|         ranging from min_y2 to max_y2                              |
|---------------<--  End of Files or Arguments Needed -->------------|
|--------------------------------------------------------------------|
|----------------<--  Start of Example and Usage -->-----------------|
| Example: %mn2_pr(alpha=.05,sides=1,y1=0,min_y2=.2,max_y2=1,        |
|                  inc_y2=.1,std=1, n1=51,n2=51,unit=gals,plot=p);   |
| Usage:   %mn2_pr(ALPHA=.05,SIDES=2,N1=.,N2=.,Y1=.,MIN_Y2=.,        |
|              MAX_Y2=.,INC_Y2=.,STD=.,PLOT= , UNIT=);               |
| Reference: Bergstralh, EJ.  SAS macros for sample size and power   |
|            calculations.  Proceedings of the 9th annual SAS Users  |
|            Group International Conference.                         |
|            Equation #8.                                            |
\-------------------<--  End of Example and Usage -->---------------*/
%MACRO  mn2_pr(ALPHA=.05,SIDES=2,N1=.,N2=.,Y1=.,MIN_Y2=.,
               MAX_Y2=.,INC_Y2=.,STD=.,PLOT= , UNIT=);
/*--------------------------------------------\
| Author:   Michael Riggs and Eric Bergstralh;|
| Created:  November 16, 1992;                |
| Modified: January 9, 1998;                  |
| Purpose:  Power for two independent means;  |
\--------------------------------------------*/
 OPTIONS MISSING=' ' NOCENTER;
   %LET PLOT=%UPCASE(&PLOT);
 DATA T1;
      ALPHA=&ALPHA;
      SIDES=&SIDES;
      Y1=&Y1;
      MIN_Y2=&MIN_Y2;
      MAX_Y2=&MAX_Y2;
      INC_Y2=&INC_Y2;
      N1=&N1;
      N2=&N2;
      STD=&STD;
    ZALPHA=(PROBIT(1-ALPHA))*(SIDES=1) + (PROBIT(1-ALPHA/2))*(SIDES=2);
    N=N1+N2;
    IF MAX_Y2=. THEN DO;
      MAX_Y2=MIN_Y2+1; INC_Y2=MIN_Y2+2;  *NEED 1 EXEC OF DO;
    END;
    TN1=N1;
    TN2=N2;
    df=tn1+tn2-2; 
    tn1=int( tn1*(df+1)/(df+3) );  *decrease n for t-test 1/9/98;    
    tn2=int( tn2*(df+1)/(df+3) );  *decrease n for t-test;    

    TY1=Y1;
    TSTD=STD;
    DO Y2=MIN_Y2 TO MAX_Y2 BY INC_Y2;
       ZBETA=ABS(Y2-TY1)/(TSTD*SQRT(1/TN1+1/TN2))-ZALPHA;
       POWER=PROBNORM(ZBETA);
       OUTPUT;
     END;
 LABEL Y1="Group 1*Mean (&unit)"
       Y2="Group 2*Mean (&unit)"
       N1='Group 1*Sample size'
       N2='Group 2*Sample size';
FOOTNOTE1 '   ';
 run;
 PROC PRINT SPLIT='*';
  ID Y1; var  Y2 POWER;

 TITLE2
   "POWER ESTIMATES FOR COMPARING TWO MEANS";
 title3
 "Alpha=&alpha, Sides=&sides, Pooled SD=&std, N1=&n1, N2=&n2";
 title4"Group 1 true mean=&y1(&unit)";
%IF "&MAX_Y2" NE "" %THEN %DO;
    %IF &PLOT= P  %THEN %DO;
      PROC PLOT NOLEGEND;
           PLOT POWER*Y2/ VAXIS=0 TO 1.0 BY .10
                          HAXIS=&MIN_Y2 TO &MAX_Y2 BY &INC_Y2;
      LABEL Y2="True Group 2 mean(&unit)";
    %END;
    %ELSE %IF &PLOT= G  %THEN %DO;
          PROC GPLOT ;
      PLOT POWER*Y2/ vaxis=0 to 1.0 by .2;
      SYMBOL1  Font= v=none i=join l=1;
      LABEL Y2="True Group 2 mean(&unit)";
   run;
    %END;
    RUN;
%END;
%MEND mn2_pr;